NewSnowAnalysis
Shihong Jia
Oct. 16, 2017
Load packages
library(ggplot2)
library(cowplot)
library(lme4)
library(pbkrtest)
library(tidyr)
library(dplyr)
library(optimx)
library(broom)Data description
Seven censuses were conducted in this experiment (from 14fa to 17fa, twice a year). Only those seedlings inside exclosures were used in the analyses. Treatment: C (no manipulation), S (snow removal treatment).
Initial analysis
getwd()## [1] "E:/R/My code/Git/Changbaishan/Changbaishan1"## Read in data
snowdat.ag_adult <- read.csv("data/snowdat.ag_adult.csv")
snowdat.ag_adult$census <- factor(snowdat.ag_adult$census)
snowdat.ag_adult$site <- factor(snowdat.ag_adult$site)
snowdat.ag_adult$census <- factor(snowdat.ag_adult$census, levels=c('15sp', '15fa', '16sp', '16fa','17sp','17fa'))
summary(snowdat.ag_adult)## X site quadrat sp. treatment
## Min. : 1.0 1:445 Min. :102 FRMA :217 C:657
## 1st Qu.: 291.2 2:455 1st Qu.:203 TIAM :189 S:505
## Median : 581.5 3:262 Median :306 PHSC :115
## Mean : 581.5 Mean :299 EUMA : 82
## 3rd Qu.: 871.8 3rd Qu.:404 DEGL : 60
## Max. :1162.0 Max. :508 ACPS : 52
## (Other):447
## growth.form census deaths total
## shrub:467 15sp:208 Min. : 0.0000 Min. : 1.000
## tree :695 15fa:117 1st Qu.: 0.0000 1st Qu.: 1.000
## 16sp: 94 Median : 0.0000 Median : 2.000
## 16fa:225 Mean : 0.6945 Mean : 2.543
## 17sp:219 3rd Qu.: 1.0000 3rd Qu.: 3.000
## 17fa:299 Max. :17.0000 Max. :21.000
##
## survs quad.unique quad.sp A.sum
## Min. : 0.000 1_1_308: 43 1_1_102_EUAL: 6 Min. :4487
## 1st Qu.: 1.000 2_1_403: 43 1_1_102_PHSC: 6 1st Qu.:5335
## Median : 1.000 2_1_501: 43 1_1_104_DEGL: 6 Median :5708
## Mean : 1.849 2_1_104: 37 1_1_104_EUMA: 6 Mean :5902
## 3rd Qu.: 2.000 1_1_507: 35 1_1_203_DEGL: 6 3rd Qu.:6533
## Max. :14.000 1_1_409: 32 1_1_205_PHSC: 6 Max. :9223
## (Other):929 (Other) :1126
## A.con A.het con.dens A.con_curt
## Min. : 0.000 Min. :2827 Min. : 0.000 Min. : 0.000
## 1st Qu.: 0.000 1st Qu.:4766 1st Qu.: 0.000 1st Qu.: 0.000
## Median : 3.047 Median :5455 Median : 1.000 Median : 1.450
## Mean : 446.019 Mean :5456 Mean : 4.998 Mean : 4.357
## 3rd Qu.: 648.325 3rd Qu.:6120 3rd Qu.: 8.000 3rd Qu.: 8.655
## Max. :4003.948 Max. :9223 Max. :41.000 Max. :15.879
##
## A.het_curt A.sum_curt con.dens_curt
## Min. :14.14 Min. :16.49 Min. :0.000
## 1st Qu.:16.83 1st Qu.:17.47 1st Qu.:0.000
## Median :17.60 Median :17.87 Median :1.000
## Mean :17.53 Mean :18.03 Mean :1.010
## 3rd Qu.:18.29 3rd Qu.:18.69 3rd Qu.:2.000
## Max. :20.97 Max. :20.97 Max. :3.448
## ###------------------------ Initial plots -------------------------###
ggplot(data=snowdat.ag_adult,
aes(x=treatment, y=survs/total, colour=as.factor(census))) +
geom_boxplot() +
facet_wrap(~growth.form)
group_by(snowdat.ag_adult, treatment, growth.form, census) %>%
summarise(surv = mean(survs/total),
se.survs=sd(survs/total)/sqrt(length(survs/total))) %>%
ggplot(aes(x=treatment, y=surv, ymin=surv-se.survs, ymax=surv+se.survs,
colour=as.factor(census))) + geom_errorbar() + geom_point() +
facet_wrap(~growth.form) 
group_by(snowdat.ag_adult, treatment, growth.form, census) %>%
summarise(surv = mean(survs/total),
se.survs=sd(survs/total)/sqrt(length(survs/total))) %>%
ggplot(aes(x=census, y=surv, ymin=surv-se.survs, ymax=surv+se.survs,
colour=as.factor(treatment))) + geom_errorbar() + geom_point() +
facet_wrap(~growth.form) 
###----------------------- Survival models ------------------------###
snowdat.ag_adult$site <- factor(snowdat.ag_adult$site)
contrasts(snowdat.ag_adult$site) <- "contr.sum"
#### tree
## with or without neighboring trees
mod.snow.tree.surv1 <- glmer(cbind(survs, deaths) ~ site + treatment + A.con_curt + A.het_curt +
(1|quad.unique) + (1|sp.) + (1|census),
data=subset(snowdat.ag_adult, growth.form=='tree'), family=binomial,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
# without neighboring trees
mod.snow.tree.surv2 <- glmer(cbind(survs, deaths) ~ site + treatment +
(1|quad.unique) + (1|sp.) + (1|census),
data=subset(snowdat.ag_adult, growth.form=='tree'), family=binomial,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
anova(mod.snow.tree.surv1, mod.snow.tree.surv2)## Data: subset(snowdat.ag_adult, growth.form == "tree")
## Models:
## mod.snow.tree.surv2: cbind(survs, deaths) ~ site + treatment + (1 | quad.unique) +
## mod.snow.tree.surv2: (1 | sp.) + (1 | census)
## mod.snow.tree.surv1: cbind(survs, deaths) ~ site + treatment + A.con_curt + A.het_curt +
## mod.snow.tree.surv1: (1 | quad.unique) + (1 | sp.) + (1 | census)
## Df AIC BIC logLik deviance Chisq Chi Df
## mod.snow.tree.surv2 7 1343.2 1375 -664.62 1329.2
## mod.snow.tree.surv1 9 1344.1 1385 -663.07 1326.1 3.0972 2
## Pr(>Chisq)
## mod.snow.tree.surv2
## mod.snow.tree.surv1 0.2125# adding area of neighboring trees make no sense to the model, but I would like to keep it as the same as the survival analysis of the pesticide experiment.
# model sensitivity tests
mod.snow.tree.surv.diags <- augment(mod.snow.tree.surv2)
qplot(data=mod.snow.tree.surv.diags, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
geom_hline(yintercept=0, linetype='dotted')## `geom_smooth()` using method = 'loess'
qplot(data=mod.snow.tree.surv.diags, x=treatment, y=.wtres, geom="boxplot") 
library(car)##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recodeqqPlot(resid(mod.snow.tree.surv1))
library(sjPlot)##
## Attaching package: 'sjPlot'## The following objects are masked from 'package:cowplot':
##
## plot_grid, save_plotsjp.lmer(mod.snow.tree.surv2, type='fe')
sjp.lmer(mod.snow.tree.surv2, type='re', sort.est="sort.all" )## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.## Plotting random effects...## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.## Plotting random effects...
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...
summary(mod.snow.tree.surv2, correlation=FALSE)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(survs, deaths) ~ site + treatment + (1 | quad.unique) +
## (1 | sp.) + (1 | census)
## Data: subset(snowdat.ag_adult, growth.form == "tree")
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",
## "Nelder-Mead")))
##
## AIC BIC logLik deviance df.resid
## 1343.2 1375.0 -664.6 1329.2 688
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.2984 -0.6223 0.2682 0.6565 3.0701
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.1182 0.3437
## sp. (Intercept) 0.7354 0.8576
## census (Intercept) 1.5785 1.2564
## Number of obs: 695, groups: quad.unique, 60; sp., 15; census, 6
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.38884 0.58816 0.661 0.509
## site1 0.14549 0.10515 1.384 0.166
## site2 -0.04758 0.10295 -0.462 0.644
## treatmentS -0.07413 0.14617 -0.507 0.612##-----------------------------##
# Snow removal treatment didn't change the survival of seedlings,
# but the BLUP values of census show that snow removal duration
# may has a potential relationship with survival.
#### duration
## duration is presented as months from the beginning of the experiment
snowdat.ag_adult$duration <- ifelse(snowdat.ag_adult$census=="15sp", 9,
ifelse(snowdat.ag_adult$census=="15fa", 12,
ifelse(snowdat.ag_adult$census=="16sp", 21,
ifelse(snowdat.ag_adult$census=="16fa", 24,
ifelse(snowdat.ag_adult$census=="17sp", 32,
ifelse(snowdat.ag_adult$census=="17fa", 35, NA))))))
summary(snowdat.ag_adult)## X site quadrat sp. treatment
## Min. : 1.0 1:445 Min. :102 FRMA :217 C:657
## 1st Qu.: 291.2 2:455 1st Qu.:203 TIAM :189 S:505
## Median : 581.5 3:262 Median :306 PHSC :115
## Mean : 581.5 Mean :299 EUMA : 82
## 3rd Qu.: 871.8 3rd Qu.:404 DEGL : 60
## Max. :1162.0 Max. :508 ACPS : 52
## (Other):447
## growth.form census deaths total
## shrub:467 15sp:208 Min. : 0.0000 Min. : 1.000
## tree :695 15fa:117 1st Qu.: 0.0000 1st Qu.: 1.000
## 16sp: 94 Median : 0.0000 Median : 2.000
## 16fa:225 Mean : 0.6945 Mean : 2.543
## 17sp:219 3rd Qu.: 1.0000 3rd Qu.: 3.000
## 17fa:299 Max. :17.0000 Max. :21.000
##
## survs quad.unique quad.sp A.sum
## Min. : 0.000 1_1_308: 43 1_1_102_EUAL: 6 Min. :4487
## 1st Qu.: 1.000 2_1_403: 43 1_1_102_PHSC: 6 1st Qu.:5335
## Median : 1.000 2_1_501: 43 1_1_104_DEGL: 6 Median :5708
## Mean : 1.849 2_1_104: 37 1_1_104_EUMA: 6 Mean :5902
## 3rd Qu.: 2.000 1_1_507: 35 1_1_203_DEGL: 6 3rd Qu.:6533
## Max. :14.000 1_1_409: 32 1_1_205_PHSC: 6 Max. :9223
## (Other):929 (Other) :1126
## A.con A.het con.dens A.con_curt
## Min. : 0.000 Min. :2827 Min. : 0.000 Min. : 0.000
## 1st Qu.: 0.000 1st Qu.:4766 1st Qu.: 0.000 1st Qu.: 0.000
## Median : 3.047 Median :5455 Median : 1.000 Median : 1.450
## Mean : 446.019 Mean :5456 Mean : 4.998 Mean : 4.357
## 3rd Qu.: 648.325 3rd Qu.:6120 3rd Qu.: 8.000 3rd Qu.: 8.655
## Max. :4003.948 Max. :9223 Max. :41.000 Max. :15.879
##
## A.het_curt A.sum_curt con.dens_curt duration
## Min. :14.14 Min. :16.49 Min. :0.000 Min. : 9.0
## 1st Qu.:16.83 1st Qu.:17.47 1st Qu.:0.000 1st Qu.:12.0
## Median :17.60 Median :17.87 Median :1.000 Median :24.0
## Mean :17.53 Mean :18.03 Mean :1.010 Mean :24.2
## 3rd Qu.:18.29 3rd Qu.:18.69 3rd Qu.:2.000 3rd Qu.:35.0
## Max. :20.97 Max. :20.97 Max. :3.448 Max. :35.0
## #### tree
mod.snow.tree.surv3 <- glmer(cbind(survs, deaths) ~ site + treatment*duration +
A.con_curt + A.het_curt +
(1|quad.unique) + (1|sp.),
data=subset(snowdat.ag_adult, growth.form=='tree'), family=binomial,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
mod.snow.tree.surv4 <- glmer(cbind(survs, deaths) ~ treatment*duration +
A.con_curt + A.het_curt +
(1|quad.unique) + (1|census) + (1|sp.),
data=subset(snowdat.ag_adult, growth.form=='tree'), family=binomial,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
anova(mod.snow.tree.surv1, mod.snow.tree.surv3, mod.snow.tree.surv4)## Data: subset(snowdat.ag_adult, growth.form == "tree")
## Models:
## mod.snow.tree.surv1: cbind(survs, deaths) ~ site + treatment + A.con_curt + A.het_curt +
## mod.snow.tree.surv1: (1 | quad.unique) + (1 | sp.) + (1 | census)
## mod.snow.tree.surv4: cbind(survs, deaths) ~ treatment * duration + A.con_curt + A.het_curt +
## mod.snow.tree.surv4: (1 | quad.unique) + (1 | census) + (1 | sp.)
## mod.snow.tree.surv3: cbind(survs, deaths) ~ site + treatment * duration + A.con_curt +
## mod.snow.tree.surv3: A.het_curt + (1 | quad.unique) + (1 | sp.)
## Df AIC BIC logLik deviance Chisq Chi Df
## mod.snow.tree.surv1 9 1344.1 1385.0 -663.07 1326.1
## mod.snow.tree.surv4 9 1285.9 1326.8 -633.93 1267.9 58.273 0
## mod.snow.tree.surv3 10 1428.2 1473.6 -704.09 1408.2 0.000 1
## Pr(>Chisq)
## mod.snow.tree.surv1
## mod.snow.tree.surv4 <2e-16 ***
## mod.snow.tree.surv3 1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod.snow.tree.surv.diags <- augment(mod.snow.tree.surv4)
qplot(data=mod.snow.tree.surv.diags, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
geom_hline(yintercept=0, linetype='dotted')## `geom_smooth()` using method = 'loess'
sjp.lmer(mod.snow.tree.surv4, type='re', sort.est="sort.all" )## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.## Plotting random effects...
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...
qqPlot(resid(mod.snow.tree.surv4))
summary(mod.snow.tree.surv4, correlation=FALSE)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## cbind(survs, deaths) ~ treatment * duration + A.con_curt + A.het_curt +
## (1 | quad.unique) + (1 | census) + (1 | sp.)
## Data: subset(snowdat.ag_adult, growth.form == "tree")
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",
## "Nelder-Mead")))
##
## AIC BIC logLik deviance df.resid
## 1285.9 1326.8 -633.9 1267.9 686
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.1175 -0.5715 0.2764 0.6530 2.6078
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.09238 0.3039
## sp. (Intercept) 0.49931 0.7066
## census (Intercept) 0.36476 0.6040
## Number of obs: 695, groups: quad.unique, 60; sp., 15; census, 6
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.10421 1.40830 -1.494 0.13514
## treatmentS -2.86182 0.43677 -6.552 5.67e-11 ***
## duration 0.08102 0.02725 2.973 0.00295 **
## A.con_curt -0.03530 0.02358 -1.497 0.13439
## A.het_curt 0.06162 0.06669 0.924 0.35552
## treatmentS:duration 0.10083 0.01482 6.804 1.02e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## only according to AIC values, model 4 is better. But the model looks problemic.Other approaches?
#### season
snowdat.ag_adult$season <- ifelse(snowdat.ag_adult$census=='15fa'|
snowdat.ag_adult$census=='16fa'|
snowdat.ag_adult$census=='17fa', 'fa', 'sp')
mod.snow.tree.surv5 <- glmer(cbind(survs, deaths) ~ site + treatment*season +
A.con_curt + A.het_curt +
(1|quad.unique) + (1|census) + (1|sp.),
data=subset(snowdat.ag_adult, growth.form=='tree'), family=binomial,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
anova(mod.snow.tree.surv1, mod.snow.tree.surv3, mod.snow.tree.surv5)## Data: subset(snowdat.ag_adult, growth.form == "tree")
## Models:
## mod.snow.tree.surv1: cbind(survs, deaths) ~ site + treatment + A.con_curt + A.het_curt +
## mod.snow.tree.surv1: (1 | quad.unique) + (1 | sp.) + (1 | census)
## mod.snow.tree.surv3: cbind(survs, deaths) ~ site + treatment * duration + A.con_curt +
## mod.snow.tree.surv3: A.het_curt + (1 | quad.unique) + (1 | sp.)
## mod.snow.tree.surv5: cbind(survs, deaths) ~ site + treatment * season + A.con_curt +
## mod.snow.tree.surv5: A.het_curt + (1 | quad.unique) + (1 | census) + (1 | sp.)
## Df AIC BIC logLik deviance Chisq Chi Df
## mod.snow.tree.surv1 9 1344.1 1385.0 -663.07 1326.1
## mod.snow.tree.surv3 10 1428.2 1473.6 -704.09 1408.2 0.00 1
## mod.snow.tree.surv5 11 1308.0 1358.0 -643.01 1286.0 122.17 1
## Pr(>Chisq)
## mod.snow.tree.surv1
## mod.snow.tree.surv3 1
## mod.snow.tree.surv5 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1sjp.lmer(mod.snow.tree.surv5, type='re', sort.est="sort.all" )## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.## Plotting random effects...
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...
summary(mod.snow.tree.surv5, correlation=FALSE)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(survs, deaths) ~ site + treatment * season + A.con_curt +
## A.het_curt + (1 | quad.unique) + (1 | census) + (1 | sp.)
## Data: subset(snowdat.ag_adult, growth.form == "tree")
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",
## "Nelder-Mead")))
##
## AIC BIC logLik deviance df.resid
## 1308 1358 -643 1286 684
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.7248 -0.6124 0.2664 0.6506 2.6048
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.09498 0.3082
## sp. (Intercept) 0.56548 0.7520
## census (Intercept) 1.03429 1.0170
## Number of obs: 695, groups: quad.unique, 60; sp., 15; census, 6
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.12850 1.47414 -0.087 0.93054
## site1 0.12602 0.11359 1.109 0.26726
## site2 -0.02943 0.10722 -0.274 0.78371
## treatmentS 0.48850 0.16922 2.887 0.00389 **
## seasonsp -0.83093 0.84608 -0.982 0.32605
## A.con_curt -0.02656 0.02628 -1.011 0.31208
## A.het_curt 0.06739 0.07128 0.946 0.34440
## treatmentS:seasonsp -1.39301 0.22668 -6.145 7.98e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Maybe the model considering season as fixed effect makes more sense?
Survival analysis
# how about the shrubs?
#### shrub
mod.snow.shrub.surv1 <- glmer(cbind(survs, deaths) ~ site + treatment + A.con_curt + A.het_curt +
(1|quad.unique) + (1|sp.) + (1|census),
data=subset(snowdat.ag_adult, growth.form=='shrub'), family=binomial,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))## fixed-effect model matrix is rank deficient so dropping 1 column / coefficientmod.snow.shrub.surv2 <- glmer(cbind(survs, deaths) ~ site + treatment*duration +
(1|quad.unique) + (1|sp.),
data=subset(snowdat.ag_adult, growth.form=='shrub'), family=binomial,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
mod.snow.shrub.surv3 <- glmer(cbind(survs, deaths) ~ site + treatment*season +
(1|quad.unique) + (1|sp.) + (1|census),
data=subset(snowdat.ag_adult, growth.form=='shrub'), family=binomial,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
mod.snow.shrub.surv4 <- glmer(cbind(survs, deaths) ~ site + treatment + season +
(1|quad.unique) + (1|sp.) + (1|census),
data=subset(snowdat.ag_adult, growth.form=='shrub'), family=binomial,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
anova(mod.snow.shrub.surv1, mod.snow.shrub.surv2, mod.snow.shrub.surv3, mod.snow.shrub.surv4)## Data: subset(snowdat.ag_adult, growth.form == "shrub")
## Models:
## mod.snow.shrub.surv1: cbind(survs, deaths) ~ site + treatment + A.con_curt + A.het_curt +
## mod.snow.shrub.surv1: (1 | quad.unique) + (1 | sp.) + (1 | census)
## mod.snow.shrub.surv2: cbind(survs, deaths) ~ site + treatment * duration + (1 | quad.unique) +
## mod.snow.shrub.surv2: (1 | sp.)
## mod.snow.shrub.surv4: cbind(survs, deaths) ~ site + treatment + season + (1 | quad.unique) +
## mod.snow.shrub.surv4: (1 | sp.) + (1 | census)
## mod.snow.shrub.surv3: cbind(survs, deaths) ~ site + treatment * season + (1 | quad.unique) +
## mod.snow.shrub.surv3: (1 | sp.) + (1 | census)
## Df AIC BIC logLik deviance Chisq Chi Df
## mod.snow.shrub.surv1 8 295.55 328.72 -139.77 279.55
## mod.snow.shrub.surv2 8 314.37 347.54 -149.19 298.37 0.0000 0
## mod.snow.shrub.surv4 8 297.13 330.30 -140.57 281.13 17.2393 0
## mod.snow.shrub.surv3 9 298.60 335.92 -140.30 280.60 0.5299 1
## Pr(>Chisq)
## mod.snow.shrub.surv1
## mod.snow.shrub.surv2 1.0000
## mod.snow.shrub.surv4 <2e-16 ***
## mod.snow.shrub.surv3 0.4666
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod.snow.shrub.surv.diags <- augment(mod.snow.shrub.surv3)
qplot(data=mod.snow.shrub.surv.diags, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
geom_hline(yintercept=0, linetype='dotted')## `geom_smooth()` using method = 'loess'
qqPlot(resid(mod.snow.shrub.surv3))
summary(mod.snow.shrub.surv3, correlation=FALSE)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## cbind(survs, deaths) ~ site + treatment * season + (1 | quad.unique) +
## (1 | sp.) + (1 | census)
## Data: subset(snowdat.ag_adult, growth.form == "shrub")
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",
## "Nelder-Mead")))
##
## AIC BIC logLik deviance df.resid
## 298.6 335.9 -140.3 280.6 458
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -9.8448 0.0851 0.1316 0.2548 1.7041
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.4779 0.6913
## sp. (Intercept) 0.6600 0.8124
## census (Intercept) 1.2459 1.1162
## Number of obs: 467, groups: quad.unique, 52; sp., 14; census, 6
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.65939 0.95875 4.860 1.17e-06 ***
## site1 0.07082 0.28253 0.251 0.8021
## site2 0.14578 0.28483 0.512 0.6088
## treatmentS -1.59249 0.75305 -2.115 0.0345 *
## seasonsp -1.14958 1.16982 -0.983 0.3258
## treatmentS:seasonsp -0.60741 0.81576 -0.745 0.4565
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##------------------------- Density dependence-------------------------##
## For all tree species together
ggplot(data=subset(snowdat.ag_adult, growth.form=='tree'),
aes(x=total, y=survs/total, colour=treatment)) +
geom_jitter()+
geom_smooth(aes(weight = total), method='glm', method.args=list(family="binomial")) 
# no density dependence for seedlings overall
snowdat.ag_adult$dens <- (snowdat.ag_adult$survs - mean(snowdat.ag_adult$survs))/
(2*sd(snowdat.ag_adult$survs))
mod.snow.tree.surv_ddst1 <- glmer(cbind(survs, deaths) ~ site + treatment*season +
A.con_curt + A.het_curt + dens +
(1|quad.unique) + (1|sp.) + (1|census),
data=subset(snowdat.ag_adult, growth.form=='tree'), family=binomial,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
mod.snow.tree.surv_ddst2 <- glmer(cbind(survs, deaths) ~ treatment*duration +
A.con_curt + A.het_curt + dens +
(1|quad.unique) + (1|sp.),
data=subset(snowdat.ag_adult, growth.form=='tree'), family=binomial,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
anova(mod.snow.tree.surv5, mod.snow.tree.surv_ddst1, mod.snow.tree.surv_ddst2)## Data: subset(snowdat.ag_adult, growth.form == "tree")
## Models:
## mod.snow.tree.surv_ddst2: cbind(survs, deaths) ~ treatment * duration + A.con_curt + A.het_curt +
## mod.snow.tree.surv_ddst2: dens + (1 | quad.unique) + (1 | sp.)
## mod.snow.tree.surv5: cbind(survs, deaths) ~ site + treatment * season + A.con_curt +
## mod.snow.tree.surv5: A.het_curt + (1 | quad.unique) + (1 | census) + (1 | sp.)
## mod.snow.tree.surv_ddst1: cbind(survs, deaths) ~ site + treatment * season + A.con_curt +
## mod.snow.tree.surv_ddst1: A.het_curt + dens + (1 | quad.unique) + (1 | sp.) + (1 |
## mod.snow.tree.surv_ddst1: census)
## Df AIC BIC logLik deviance Chisq Chi Df
## mod.snow.tree.surv_ddst2 9 1263.2 1304.1 -622.61 1245.2
## mod.snow.tree.surv5 11 1308.0 1358.0 -643.01 1286.0 0.000 2
## mod.snow.tree.surv_ddst1 12 1229.2 1283.7 -602.59 1205.2 80.829 1
## Pr(>Chisq)
## mod.snow.tree.surv_ddst2
## mod.snow.tree.surv5 1
## mod.snow.tree.surv_ddst1 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod.snow.treeddst.surv.diags <- augment(mod.snow.tree.surv_ddst1)
qplot(data=mod.snow.treeddst.surv.diags, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
geom_hline(yintercept=0, linetype='dotted')## `geom_smooth()` using method = 'loess'
qqPlot(resid(mod.snow.tree.surv_ddst1))
summary(mod.snow.tree.surv_ddst1, correlation=FALSE)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(survs, deaths) ~ site + treatment * season + A.con_curt +
## A.het_curt + dens + (1 | quad.unique) + (1 | sp.) + (1 | census)
## Data: subset(snowdat.ag_adult, growth.form == "tree")
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",
## "Nelder-Mead")))
##
## AIC BIC logLik deviance df.resid
## 1229.2 1283.7 -602.6 1205.2 683
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.3732 -0.5405 0.3206 0.6915 2.4750
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.08275 0.2877
## sp. (Intercept) 0.53373 0.7306
## census (Intercept) 0.58770 0.7666
## Number of obs: 695, groups: quad.unique, 60; sp., 15; census, 6
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.80072 1.41077 -0.568 0.57032
## site1 0.21745 0.11167 1.947 0.05151 .
## site2 -0.17298 0.10577 -1.635 0.10196
## treatmentS 0.52767 0.16768 3.147 0.00165 **
## seasonsp -0.50795 0.64769 -0.784 0.43290
## A.con_curt -0.02498 0.02609 -0.957 0.33840
## A.het_curt 0.10022 0.07068 1.418 0.15623
## dens 0.88569 0.10374 8.537 < 2e-16 ***
## treatmentS:seasonsp -1.24057 0.23146 -5.360 8.33e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##------------------------- individual species -----------------------------##
## And now pulling the species apart.
## First need to filter down to species that can be analysed - need enough
## reps (lets say 10) and variation in abundance (at least >5
sp.snow.surv.counts <- summarise(group_by(snowdat.ag_adult, sp.),
n.quadrat=length(survs), max.dens=max(survs), min.dens=min(survs))
sp.snow.surv.counts$range <- sp.snow.surv.counts$max.dens - sp.snow.surv.counts$min.dens
sp.snow.surv.sel <- sp.snow.surv.counts$sp.[sp.snow.surv.counts$n.quadrat > 10 & sp.snow.surv.counts$range >5]
sp.snow.surv.counts[sp.snow.surv.counts$sp. %in% sp.snow.surv.sel,]## # A tibble: 5 × 5
## sp. n.quadrat max.dens min.dens range
## <fctr> <int> <int> <int> <int>
## 1 ABNE 17 9 0 9
## 2 ACPS 52 6 0 6
## 3 EUVE 39 6 0 6
## 4 FRMA 217 14 0 14
## 5 TIAM 189 11 0 11ggplot(data=subset(snowdat.ag_adult, sp. %in% sp.snow.surv.sel),
aes(x=total, y=survs/total, colour=treatment)) +
geom_jitter()+
geom_smooth(aes(weight = total), method='glm',
method.args=list(family="binomial"), fullrange=T) +
facet_wrap(~ sp., scales='free')## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
# it seems only six species have enough data, but conference intervals of EUVE are too huge,
# ABNE and ACTE are excluded becasue meeting issues when fitting models
## quick plot of only the species for which there appear to be enough data
sp.snow.surv.sel2 <- c('ACPS','FRMA','TIAM')
snow.survmods <- sapply(sp.snow.surv.sel2, function(sp){
print(sp)
spdat <- filter(snowdat.ag_adult, sp. == sp)
spdat <- droplevels(spdat)
mod <- glmer(cbind(survs, deaths) ~ site + treatment*season + dens +
A.con_curt + A.het_curt + (1|quad.unique) + (1|census),
data=spdat, family=binomial,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead')))
)
}, simplify=FALSE)## [1] "ACPS"## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : unable to evaluate scaled gradient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge: degenerate Hessian with 1 negative
## eigenvalues## [1] "FRMA"
## [1] "TIAM"## diagnostic plots
indmods.s.surv.resids <- lapply(snow.survmods, augment)
indmods.s.surv.resids <- do.call('rbind', mapply(function(spnm, dat) {
dat$species <- spnm
return(dat)
}, spnm=as.list(names(indmods.s.surv.resids)), dat=indmods.s.surv.resids,
SIMPLIFY=FALSE))
ggplot(data=indmods.s.surv.resids, aes(x=.fitted, y=.wtres, colour=treatment)) +
geom_jitter() + geom_smooth() + facet_wrap(~species, scales='free') +
geom_hline(yintercept=0, linetype='dashed')## `geom_smooth()` using method = 'loess'
library(lattice)
lapply(snow.survmods, function(mod){
dotplot(ranef(mod, condVar=TRUE))})## $ACPS
## $ACPS$quad.unique
##
## $ACPS$census
##
##
## $FRMA
## $FRMA$quad.unique
##
## $FRMA$census
##
##
## $TIAM
## $TIAM$quad.unique
##
## $TIAM$census
lapply(snow.survmods, function(x) summary(x, correlation=FALSE))## Warning in vcov.merMod(object, use.hessian = use.hessian): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Warning in vcov.merMod(object, correlation = correlation, sigm = sig): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## $ACPS
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## cbind(survs, deaths) ~ site + treatment * season + dens + A.con_curt +
## A.het_curt + (1 | quad.unique) + (1 | census)
## Data: spdat
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",
## "Nelder-Mead")))
##
## AIC BIC logLik deviance df.resid
## 75.5 97.0 -26.7 53.5 41
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.5213 -0.3807 0.1999 0.5332 2.0530
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 8.874e-05 0.00942
## census (Intercept) 1.755e+00 1.32469
## Number of obs: 52, groups: quad.unique, 26; census, 6
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 19.1402 8.9994 2.127 0.0334 *
## site2 -1.2808 0.7986 -1.604 0.1088
## site3 0.9032 1.3712 0.659 0.5101
## treatmentS 19.8869 8325.8311 0.002 0.9981
## seasonsp -2.4824 1.5969 -1.554 0.1201
## dens 0.9190 0.8284 1.109 0.2672
## A.con_curt -0.2111 0.6469 -0.326 0.7442
## A.het_curt -0.8704 0.4864 -1.790 0.0735 .
## treatmentS:seasonsp -19.5511 8325.8311 -0.002 0.9981
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## unable to evaluate scaled gradient
## Model failed to converge: degenerate Hessian with 1 negative eigenvalues
##
##
## $FRMA
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## cbind(survs, deaths) ~ site + treatment * season + dens + A.con_curt +
## A.het_curt + (1 | quad.unique) + (1 | census)
## Data: spdat
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",
## "Nelder-Mead")))
##
## AIC BIC logLik deviance df.resid
## 488.5 525.7 -233.3 466.5 206
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.7912 -0.6794 0.3483 0.7508 2.2526
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.0878 0.2963
## census (Intercept) 0.3931 0.6270
## Number of obs: 217, groups: quad.unique, 54; census, 6
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.48326 1.86984 -0.793 0.42763
## site2 -1.24517 0.38887 -3.202 0.00136 **
## site3 -0.70858 0.32464 -2.183 0.02906 *
## treatmentS 0.51463 0.26828 1.918 0.05508 .
## seasonsp -0.20714 0.56601 -0.366 0.71439
## dens 0.90654 0.16632 5.451 5.02e-08 ***
## A.con_curt 0.05652 0.03820 1.479 0.13905
## A.het_curt 0.14065 0.10067 1.397 0.16237
## treatmentS:seasonsp -1.06678 0.35881 -2.973 0.00295 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $TIAM
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## cbind(survs, deaths) ~ site + treatment * season + dens + A.con_curt +
## A.het_curt + (1 | quad.unique) + (1 | census)
## Data: spdat
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",
## "Nelder-Mead")))
##
## AIC BIC logLik deviance df.resid
## 408.3 444.0 -193.2 386.3 178
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.70668 -0.58270 -0.06585 0.76267 2.07260
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.07743 0.2783
## census (Intercept) 0.59358 0.7704
## Number of obs: 189, groups: quad.unique, 59; census, 6
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.446536 2.134263 0.678 0.498
## site2 0.130718 0.257250 0.508 0.611
## site3 -0.079597 0.246263 -0.323 0.747
## treatmentS 0.110885 0.231615 0.479 0.632
## seasonsp -1.053821 0.806889 -1.306 0.192
## dens 1.192432 0.206377 5.778 7.56e-09 ***
## A.con_curt -0.140947 0.078349 -1.799 0.072 .
## A.het_curt 0.001488 0.102478 0.015 0.988
## treatmentS:seasonsp -0.456131 0.378106 -1.206 0.228
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lapply(snow.survmods, function(x) anova(x))## $ACPS
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 2.7392 1.3696 1.3696
## treatment 1 0.7594 0.7594 0.7594
## season 1 2.8410 2.8410 2.8410
## dens 1 0.6961 0.6961 0.6961
## A.con_curt 1 0.2729 0.2729 0.2729
## A.het_curt 1 3.2023 3.2023 3.2023
## treatment:season 1 0.0000 0.0000 0.0000
##
## $FRMA
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 2.750 1.375 1.3749
## treatment 1 0.705 0.705 0.7046
## season 1 3.322 3.322 3.3221
## dens 1 32.034 32.034 32.0340
## A.con_curt 1 0.797 0.797 0.7965
## A.het_curt 1 1.471 1.471 1.4707
## treatment:season 1 8.849 8.849 8.8485
##
## $TIAM
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 0.348 0.174 0.1740
## treatment 1 0.034 0.034 0.0342
## season 1 5.277 5.277 5.2770
## dens 1 41.220 41.220 41.2204
## A.con_curt 1 3.807 3.807 3.8070
## A.het_curt 1 0.001 0.001 0.0005
## treatment:season 1 1.471 1.471 1.4712##---------------- Survival of other tree seedlings-------------------##
# sp.sel2 <- c('ACBA', 'ACPS', 'ACTE', 'FRMA', 'PIKO', 'TIAM')
## other tree species
snowdat.ag_adult.t <- subset(snowdat.ag_adult, growth.form == 'tree')
snowdat.ag_adult.t1 <- snowdat.ag_adult.t[!(snowdat.ag_adult.t$sp. %in% sp.snow.surv.sel2),]
summary(snowdat.ag_adult.t1)## X site quadrat sp. treatment
## Min. : 5.0 1:112 Min. :102.0 ACTE :51 C:121
## 1st Qu.: 144.0 2: 97 1st Qu.:203.0 ACBA :45 S:116
## Median : 303.0 3: 28 Median :306.0 ACMO :38
## Mean : 349.4 Mean :292.5 ACMA :28
## 3rd Qu.: 499.0 3rd Qu.:403.0 PIKO :23
## Max. :1158.0 Max. :508.0 ABNE :17
## (Other):35
## growth.form census deaths total survs
## shrub: 0 15sp:64 Min. :0.0000 Min. : 1.00 Min. :0.0000
## tree :237 15fa:18 1st Qu.:0.0000 1st Qu.: 1.00 1st Qu.:0.0000
## 16sp:12 Median :0.0000 Median : 1.00 Median :1.0000
## 16fa:38 Mean :0.5781 Mean : 1.57 Mean :0.9916
## 17sp:33 3rd Qu.:1.0000 3rd Qu.: 2.00 3rd Qu.:1.0000
## 17fa:72 Max. :7.0000 Max. :12.00 Max. :9.0000
##
## quad.unique quad.sp A.sum A.con
## 1_1_308: 18 1_1_305_ACMA: 6 Min. :4487 Min. : 0.000
## 1_1_109: 12 1_1_308_ACBA: 6 1st Qu.:5335 1st Qu.: 1.712
## 1_1_205: 10 1_1_308_ACMO: 6 Median :5657 Median : 70.242
## 1_1_305: 10 1_1_308_ACTE: 6 Mean :5811 Mean : 292.147
## 2_1_206: 10 2_1_102_ACMO: 6 3rd Qu.:6261 3rd Qu.: 315.070
## 3_1_109: 10 2_1_104_ACTE: 6 Max. :7822 Max. :2670.429
## (Other):167 (Other) :201
## A.het con.dens A.con_curt A.het_curt
## Min. :2906 Min. : 0.000 Min. : 0.000 Min. :14.27
## 1st Qu.:5085 1st Qu.: 1.000 1st Qu.: 1.196 1st Qu.:17.20
## Median :5450 Median : 4.000 Median : 4.126 Median :17.60
## Mean :5519 Mean : 5.738 Mean : 4.468 Mean :17.62
## 3rd Qu.:6122 3rd Qu.: 8.000 3rd Qu.: 6.805 3rd Qu.:18.29
## Max. :7779 Max. :36.000 Max. :13.874 Max. :19.81
##
## A.sum_curt con.dens_curt duration season
## Min. :16.49 Min. :0.000 Min. : 9.00 Length:237
## 1st Qu.:17.47 1st Qu.:1.000 1st Qu.: 9.00 Class :character
## Median :17.82 Median :1.587 Median :24.00 Mode :character
## Mean :17.94 Mean :1.446 Mean :23.34
## 3rd Qu.:18.43 3rd Qu.:2.000 3rd Qu.:35.00
## Max. :19.85 Max. :3.302 Max. :35.00
##
## dens
## Min. :-0.4615
## 1st Qu.:-0.4615
## Median :-0.2118
## Mean :-0.2139
## 3rd Qu.:-0.2118
## Max. : 1.7853
## mod.snow.resttree.surv <- glmer(cbind(survs, deaths) ~ site + treatment*season +
A.con_curt + A.het_curt+
(1|quad.unique) + (1|census) + (1|sp.),
data=snowdat.ag_adult.t1, family=binomial,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))## Warning in optimx.check(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower, : Parameters or bounds appear to have different scalings.
## This can cause poor performance in optimization.
## It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.mod.snow.resttree.surv.diags <- augment(mod.snow.resttree.surv)
qplot(data=mod.snow.resttree.surv.diags, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
geom_hline(yintercept=0, linetype='dotted')## `geom_smooth()` using method = 'loess'
qplot(data=mod.snow.resttree.surv.diags, x=treatment, y=.wtres, geom="boxplot") 
## plot fixed and random effects
sjp.lmer(mod.snow.resttree.surv, type='fe', p.kr = FALSE)
sjp.lmer(mod.snow.resttree.surv, type='re', sort.est="sort.all", p.kr = FALSE)## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.## Plotting random effects...## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.## Plotting random effects...
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...
qqPlot(resid(mod.snow.resttree.surv))
summary(mod.snow.resttree.surv, correlation=FALSE)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(survs, deaths) ~ site + treatment * season + A.con_curt +
## A.het_curt + (1 | quad.unique) + (1 | census) + (1 | sp.)
## Data: snowdat.ag_adult.t1
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",
## "Nelder-Mead")))
##
## AIC BIC logLik deviance df.resid
## 282.8 321.0 -130.4 260.8 226
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.8257 -0.2492 0.1888 0.4405 3.4835
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.2952 0.5433
## sp. (Intercept) 0.8839 0.9402
## census (Intercept) 1.3861 1.1773
## Number of obs: 237, groups: quad.unique, 48; sp., 12; census, 6
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.9312000 4.0089492 -0.981 0.32679
## site1 0.3240740 0.3278558 0.988 0.32292
## site2 -0.5856986 0.3843514 -1.524 0.12754
## treatmentS 1.7165958 0.5280091 3.251 0.00115 **
## seasonsp -0.0555158 1.0901599 -0.051 0.95939
## A.con_curt -0.0005923 0.0979687 -0.006 0.99518
## A.het_curt 0.2795559 0.2176052 1.285 0.19890
## treatmentS:seasonsp -4.4869496 0.8002280 -5.607 2.06e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Parameters or bounds appear to have different scalings.
## This can cause poor performance in optimization.
## It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.Recruitment analysis
###-------------------- recruitment annlysis -----------------------------###
snow.rect.ag_adult <- read.csv("data/snow.rect.ag_adult.csv")
snow.rect.ag_adult$site <- factor(snow.rect.ag_adult$site)
contrasts(snow.rect.ag_adult$site) <- "contr.sum"
## Trees
snow.t.rect_adult <- snow.rect.ag_adult[snow.rect.ag_adult$growth.form == 'tree',]
mod.snow.recr1 <- glmer(recruits ~ site + treatment*census + A.con_curt + A.het_curt +
(1|quad.unique) + (1|sp.),
data=snow.t.rect_adult, family=poisson,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
mod.snow.recr2 <- glmer(recruits ~ site + treatment + A.con_curt + A.het_curt +
(1|quad.unique) + (1|sp.),
data=snow.t.rect_adult, family=poisson,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
anova(mod.snow.recr1, mod.snow.recr2)## Data: snow.t.rect_adult
## Models:
## mod.snow.recr2: recruits ~ site + treatment + A.con_curt + A.het_curt + (1 |
## mod.snow.recr2: quad.unique) + (1 | sp.)
## mod.snow.recr1: recruits ~ site + treatment * census + A.con_curt + A.het_curt +
## mod.snow.recr1: (1 | quad.unique) + (1 | sp.)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mod.snow.recr2 8 1426.1 1456.2 -705.07 1410.1
## mod.snow.recr1 12 1348.3 1393.3 -662.16 1324.3 85.825 4 < 2.2e-16
##
## mod.snow.recr2
## mod.snow.recr1 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1sjp.lmer(mod.snow.recr1, type='re', sort.est="sort.all" )## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.## Plotting random effects...## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.## Plotting random effects...
## Random effects ofABNE, FRMA and TIAM go extramely, try to put group them against other species as the fixed effect
snow.t.rect_adult$sp.group <- ifelse(snow.t.rect_adult$sp.=='ABNE' |
snow.t.rect_adult$sp.=='FRMA' | snow.t.rect_adult$sp.=='TIAM', 1, 0)
snow.t.rect_adult$sp.group <- as.factor(snow.t.rect_adult$sp.group)
mod.snow.recr3 <- glmer(recruits ~ site + treatment*census +
A.con_curt + A.het_curt + sp.group +
(1|quad.unique) + (1|sp.),
data=snow.t.rect_adult, family=poisson,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
anova(mod.snow.recr1, mod.snow.recr2, mod.snow.recr3)## Data: snow.t.rect_adult
## Models:
## mod.snow.recr2: recruits ~ site + treatment + A.con_curt + A.het_curt + (1 |
## mod.snow.recr2: quad.unique) + (1 | sp.)
## mod.snow.recr1: recruits ~ site + treatment * census + A.con_curt + A.het_curt +
## mod.snow.recr1: (1 | quad.unique) + (1 | sp.)
## mod.snow.recr3: recruits ~ site + treatment * census + A.con_curt + A.het_curt +
## mod.snow.recr3: sp.group + (1 | quad.unique) + (1 | sp.)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mod.snow.recr2 8 1426.1 1456.2 -705.07 1410.1
## mod.snow.recr1 12 1348.3 1393.3 -662.16 1324.3 85.825 4 < 2.2e-16
## mod.snow.recr3 13 1333.1 1381.9 -653.57 1307.1 17.174 1 3.412e-05
##
## mod.snow.recr2
## mod.snow.recr1 ***
## mod.snow.recr3 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod.snow.recr.diags <- augment(mod.snow.recr1)
qplot(data=mod.snow.recr.diags, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
geom_hline(yintercept=0, linetype='dotted')## `geom_smooth()` using method = 'loess'
qplot(data = mod.snow.recr.diags, x=treatment, y=.wtres, geom="boxplot") 
sjp.lmer(mod.snow.recr3, type='fe')
sjp.lmer(mod.snow.recr3, type='re', sort.est="sort.all" )## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.## Plotting random effects...
qqPlot(resid(mod.snow.recr3))
summary(mod.snow.recr3, correlation=FALSE)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: recruits ~ site + treatment * census + A.con_curt + A.het_curt +
## sp.group + (1 | quad.unique) + (1 | sp.)
## Data: snow.t.rect_adult
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",
## "Nelder-Mead")))
##
## AIC BIC logLik deviance df.resid
## 1333.1 1381.9 -653.6 1307.1 302
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9358 -0.7704 -0.2309 0.4938 5.0218
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.05735 0.2395
## sp. (Intercept) 0.03437 0.1854
## Number of obs: 315, groups: quad.unique, 60; sp., 14
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.986064 0.866612 1.138 0.25519
## site1 -0.076718 0.070908 -1.082 0.27927
## site2 0.118672 0.065411 1.814 0.06964 .
## treatmentS -1.014233 0.329974 -3.074 0.00211 **
## census16sp 0.640176 0.150175 4.263 2.02e-05 ***
## census17sp 0.183148 0.164943 1.110 0.26684
## A.con_curt 0.008625 0.014882 0.580 0.56220
## A.het_curt -0.062159 0.046168 -1.346 0.17818
## sp.group1 1.018866 0.172732 5.899 3.67e-09 ***
## treatmentS:census16sp 0.911925 0.333331 2.736 0.00622 **
## treatmentS:census17sp 0.916593 0.343051 2.672 0.00754 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1###---------------- individual species analysis ----------------------------###
## reps (lets say 10) and variation in abundance (at least >5
sp.snowrect.counts <- summarise(group_by(snow.t.rect_adult, sp.),
n.quadrat=length(recruits), max.rect=max(recruits), min.rect=min(recruits))
sp.snowrect.counts$range <- sp.snowrect.counts$max.rect - sp.snowrect.counts$min.rect
sp.snowrect.sel <- sp.snowrect.counts$sp.[sp.snowrect.counts$n.quadrat > 10 & sp.snowrect.counts$range >5]
sp.snowrect.counts[sp.snowrect.counts$sp. %in% sp.snowrect.sel,]## # A tibble: 4 × 5
## sp. n.quadrat max.rect min.rect range
## <fctr> <int> <int> <int> <int>
## 1 ABNE 16 12 1 11
## 2 ACPS 30 9 1 8
## 3 FRMA 99 19 1 18
## 4 TIAM 104 12 1 11sp.snowrect.sel2 <- c('ACPS','FRMA','TIAM')
snowrectmods <- sapply(sp.snowrect.sel2, function(sp){
print(sp)
spdat <- filter(snow.t.rect_adult, sp. == sp)
spdat <- droplevels(spdat)
mod <- glmer(recruits ~ site + treatment*census +
A.con_curt + A.het_curt + (1|quad.unique),
data=spdat, family=poisson,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead')))
)
}, simplify=FALSE)## [1] "ACPS"## fixed-effect model matrix is rank deficient so dropping 2 columns / coefficients## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : unable to evaluate scaled gradient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge: degenerate Hessian with 1 negative
## eigenvalues## [1] "FRMA"
## [1] "TIAM"## diagnostic plots
indmods.snowrect.resids <- lapply(snowrectmods, augment)
indmods.snowrect.resids <- do.call('rbind', mapply(function(spnm, dat) {
dat$species <- spnm
return(dat)
}, spnm=as.list(names(indmods.snowrect.resids)), dat=indmods.snowrect.resids,
SIMPLIFY=FALSE))
ggplot(data=indmods.snowrect.resids, aes(x=.fitted, y=.wtres, colour=treatment)) +
geom_jitter() + geom_smooth() + facet_wrap(~species, scales='free') +
geom_hline(yintercept=0, linetype='dashed')## `geom_smooth()` using method = 'loess'## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.0033838## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.059818## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 0.43823## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small.
## fewer data values than degrees of freedom.## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at -0.0033838## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 0.059818## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 0## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 0.43823
lapply(snowrectmods, function(mod){
dotplot(ranef(mod, condVar=TRUE))})## $ACPS
## $ACPS$quad.unique
##
##
## $FRMA
## $FRMA$quad.unique
##
##
## $TIAM
## $TIAM$quad.unique
lapply(snowrectmods, function(x) summary(x, correlation=FALSE))## $ACPS
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: recruits ~ site + treatment * census + A.con_curt + A.het_curt +
## (1 | quad.unique)
## Data: spdat
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",
## "Nelder-Mead")))
##
## AIC BIC logLik deviance df.resid
## 112.2 124.8 -47.1 94.2 21
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.3579 -0.5173 -0.1187 0.3802 2.6937
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.0003189 0.01786
## Number of obs: 30, groups: quad.unique, 24
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.5180 3.7233 -0.945 0.345
## site2 0.3474 0.3244 1.071 0.284
## site3 -0.3199 0.5520 -0.580 0.562
## treatmentS -0.1368 0.6354 -0.215 0.830
## census16sp 0.4643 1.1322 0.410 0.682
## census17sp 1.0981 1.1792 0.931 0.352
## A.con_curt 0.2011 0.2257 0.891 0.373
## A.het_curt 0.1149 0.1944 0.591 0.554
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 2 columns / coefficients
## convergence code: 0
## unable to evaluate scaled gradient
## Model failed to converge: degenerate Hessian with 1 negative eigenvalues
##
##
## $FRMA
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: recruits ~ site + treatment * census + A.con_curt + A.het_curt +
## (1 | quad.unique)
## Data: spdat
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",
## "Nelder-Mead")))
##
## AIC BIC logLik deviance df.resid
## 438.8 467.3 -208.4 416.8 88
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7627 -0.8554 -0.1045 0.6032 3.0342
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.02111 0.1453
## Number of obs: 99, groups: quad.unique, 54
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.08221 1.03230 2.017 0.043690 *
## site2 0.99959 0.18393 5.435 5.49e-08 ***
## site3 0.45709 0.16984 2.691 0.007118 **
## treatmentS -1.13180 0.38565 -2.935 0.003338 **
## census16sp 0.62764 0.16391 3.829 0.000129 ***
## census17sp -0.73662 0.33055 -2.228 0.025848 *
## A.con_curt -0.06822 0.01976 -3.452 0.000556 ***
## A.het_curt -0.06138 0.05608 -1.094 0.273769
## treatmentS:census16sp 1.09505 0.40032 2.735 0.006229 **
## treatmentS:census17sp 1.28666 0.51631 2.492 0.012701 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $TIAM
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: recruits ~ site + treatment * census + A.con_curt + A.het_curt +
## (1 | quad.unique)
## Data: spdat
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",
## "Nelder-Mead")))
##
## AIC BIC logLik deviance df.resid
## 455.3 484.4 -216.7 433.3 93
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.2675 -0.5815 -0.1436 0.3405 2.2702
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.1564 0.3954
## Number of obs: 104, groups: quad.unique, 58
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.920613 1.667807 -0.552 0.58096
## site2 -0.103781 0.194799 -0.533 0.59420
## site3 -0.017650 0.194939 -0.090 0.92786
## treatmentS 0.041453 0.905515 0.046 0.96349
## census16sp 1.509932 0.529904 2.849 0.00438 **
## census17sp 0.852936 0.538295 1.584 0.11308
## A.con_curt 0.087065 0.061449 1.417 0.15652
## A.het_curt -0.002372 0.078127 -0.030 0.97578
## treatmentS:census16sp -0.175931 0.910975 -0.193 0.84686
## treatmentS:census17sp 0.267613 0.919976 0.291 0.77113
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lapply(snowrectmods, function(x) anova(x))## $ACPS
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 1.47123 0.73561 0.7356
## treatment 1 1.97066 1.97066 1.9707
## census 2 2.85826 1.42913 1.4291
## A.con_curt 1 0.89657 0.89657 0.8966
## A.het_curt 1 0.34909 0.34909 0.3491
##
## $FRMA
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 13.021 6.510 6.5105
## treatment 1 0.741 0.741 0.7409
## census 2 72.516 36.258 36.2582
## A.con_curt 1 10.690 10.690 10.6902
## A.het_curt 1 1.313 1.313 1.3128
## treatment:census 2 8.029 4.014 4.0144
##
## $TIAM
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 1.3924 0.6962 0.6962
## treatment 1 0.0314 0.0314 0.0314
## census 2 22.5802 11.2901 11.2901
## A.con_curt 1 2.3351 2.3351 2.3351
## A.het_curt 1 0.0000 0.0000 0.0000
## treatment:census 2 3.9523 1.9762 1.9762###----------------- Pooling other species -------------------------###
snow.t.rect.other_adult <- snow.t.rect_adult[!(snow.t.rect_adult$sp. %in% sp.snowrect.sel2),]
summary(snow.t.rect.other_adult)## X site quadrat sp. treatment
## Min. : 3.00 1:32 Min. :102.0 ABNE :16 C:30
## 1st Qu.: 81.25 2:37 1st Qu.:203.0 ACBA :15 S:52
## Median :142.00 3:13 Median :305.0 ACMO :12
## Mean :153.76 Mean :290.1 QUMO :10
## 3rd Qu.:225.75 3rd Qu.:404.0 ACTE : 8
## Max. :343.00 Max. :508.0 PIKO : 7
## (Other):14
## growth.form census recruits quad.unique quad.sp
## shrub: 0 15sp: 2 Min. : 1.000 1_1_205: 5 1_1_205_ACBA: 2
## tree :82 16sp:30 1st Qu.: 1.000 1_1_508: 5 1_1_409_ACBA: 2
## 17sp:50 Median : 1.000 1_1_109: 4 1_1_508_ACBA: 2
## Mean : 1.829 1_1_502: 4 2_1_203_ABNE: 2
## 3rd Qu.: 2.000 2_1_203: 4 2_1_501_ACTE: 2
## Max. :12.000 3_1_109: 4 3_1_109_ACBA: 2
## (Other):56 (Other) :70
## A.sum A.con A.het con.dens
## Min. :4487 Min. : 0.00 Min. :2906 Min. : 0.000
## 1st Qu.:5154 1st Qu.: 15.69 1st Qu.:4530 1st Qu.: 1.000
## Median :5627 Median : 191.49 Median :5267 Median : 5.500
## Mean :5699 Mean : 427.85 Mean :5271 Mean : 8.451
## 3rd Qu.:6261 3rd Qu.: 353.46 3rd Qu.:5955 3rd Qu.: 9.750
## Max. :7822 Max. :2670.43 Max. :7465 Max. :36.000
##
## A.con_curt A.het_curt A.sum_curt con.dens_curt sp.group
## Min. : 0.000 Min. :14.27 Min. :16.49 Min. :0.000 0:66
## 1st Qu.: 2.489 1st Qu.:16.55 1st Qu.:17.27 1st Qu.:1.000 1:16
## Median : 5.764 Median :17.40 Median :17.79 Median :1.764
## Mean : 5.538 Mean :17.34 Mean :17.83 Mean :1.672
## 3rd Qu.: 7.070 3rd Qu.:18.13 3rd Qu.:18.43 3rd Qu.:2.136
## Max. :13.874 Max. :19.54 Max. :19.85 Max. :3.302
## mod.snow.recr.other <- glmer(recruits ~ site + treatment*census + A.con_curt + A.het_curt +
(1|quad.unique) + (1|sp.),
data=snow.t.rect.other_adult, family=poisson,
control=glmerControl(optimizer="optimx",
optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
mod.snow.recr.other.diags <- augment(mod.snow.recr.other)
qplot(data=mod.snow.recr.other.diags, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
geom_hline(yintercept=0, linetype='dotted')## `geom_smooth()` using method = 'loess'
qplot(data = mod.snow.recr.other.diags, x=treatment, y=.wtres, geom="boxplot") 
qqPlot(resid(mod.snow.recr.other))
sjp.lmer(mod.snow.recr.other, type='fe')
sjp.lmer(mod.snow.recr.other, type='re', sort.est="sort.all" )## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.## Plotting random effects...
summary(mod.snow.recr.other, correlation=FALSE)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: recruits ~ site + treatment * census + A.con_curt + A.het_curt +
## (1 | quad.unique) + (1 | sp.)
## Data: snow.t.rect.other_adult
## Control:
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",
## "Nelder-Mead")))
##
## AIC BIC logLik deviance df.resid
## 268.3 297.2 -122.2 244.3 70
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9176 -0.3487 -0.1552 0.1270 3.4479
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.005358 0.0732
## sp. (Intercept) 0.157223 0.3965
## Number of obs: 82, groups: quad.unique, 42; sp., 11
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.93467 1.90000 0.492 0.623
## site1 0.08930 0.16227 0.550 0.582
## site2 0.04132 0.20385 0.203 0.839
## treatmentS 0.11545 1.41667 0.082 0.935
## census16sp 1.08750 1.06984 1.016 0.309
## census17sp 1.62865 1.01503 1.604 0.109
## A.con_curt 0.00336 0.04263 0.079 0.937
## A.het_curt -0.10866 0.08974 -1.211 0.226
## treatmentS:census16sp -0.16913 1.47120 -0.115 0.908
## treatmentS:census17sp -0.44983 1.43142 -0.314 0.753Growth analysis
###----------------------- Growth models ------------------------------###
snowdat.g_adult <- read.csv("data/snowdat.g_adult.csv")
snowdat.g_adult$site <- factor(snowdat.g_adult$site)
contrasts(snowdat.g_adult$site) <- "contr.sum"
#### duration
## duration is presented as months from the beginning of the experiment
snowdat.g_adult$duration <- ifelse(snowdat.g_adult$census=="15sp", 9,
ifelse(snowdat.g_adult$census=="15fa", 12,
ifelse(snowdat.g_adult$census=="16sp", 21,
ifelse(snowdat.g_adult$census=="16fa", 24,
ifelse(snowdat.g_adult$census=="17sp", 32,
ifelse(snowdat.g_adult$census=="17fa", 35, NA))))))
## season
snowdat.g_adult$season <- ifelse(snowdat.g_adult$census=='15fa'|
snowdat.g_adult$census=='16fa'|
snowdat.g_adult$census=='17fa', 'fa', 'sp')
summary(snowdat.g_adult)## X sp. site quadrat tag
## Min. : 81 FRMA :793 1:695 Min. :102.0 Min. : 1001
## 1st Qu.:5138 TIAM :379 2:926 1st Qu.:203.0 1st Qu.: 3008
## Median :6489 PHSC :179 3:462 Median :306.0 Median : 5031
## Mean :6114 EUMA : 96 Mean :295.9 Mean : 5513
## 3rd Qu.:8017 ACPS : 75 3rd Qu.:404.0 3rd Qu.: 8020
## Max. :9210 EUVE : 69 Max. :508.0 Max. :10051
## (Other):492
## treatment growth.form census grow.snd quad.unique
## C:1211 shrub: 612 15fa:142 Min. :-2.48608 2_1_210: 88
## S: 872 tree :1471 15sp:120 1st Qu.:-0.42166 2_1_403: 83
## 16fa:589 Median : 0.08960 2_1_501: 76
## 16sp:117 Mean : 0.03001 1_1_409: 69
## 17fa:693 3rd Qu.: 0.44315 2_1_104: 68
## 17sp:422 Max. : 2.53909 2_1_401: 67
## (Other):1632
## quad.sp A.sum A.con A.het
## 2_1_210_FRMA: 50 Min. :4487 Min. : 0.0 Min. :2827
## 2_1_102_FRMA: 38 1st Qu.:5335 1st Qu.: 0.0 1st Qu.:4573
## 2_1_403_FRMA: 35 Median :5708 Median : 258.6 Median :5373
## 2_1_507_FRMA: 34 Mean :5899 Mean : 577.9 Mean :5321
## 2_1_308_FRMA: 33 3rd Qu.:6533 3rd Qu.: 981.2 3rd Qu.:5974
## 2_1_401_FRMA: 28 Max. :9223 Max. :4003.9 Max. :9223
## (Other) :1865
## con.dens A.con_curt A.het_curt A.sum_curt
## Min. : 0.000 Min. : 0.000 Min. :14.14 Min. :16.49
## 1st Qu.: 0.000 1st Qu.: 0.000 1st Qu.:16.60 1st Qu.:17.47
## Median : 2.000 Median : 6.371 Median :17.51 Median :17.87
## Mean : 5.039 Mean : 5.469 Mean :17.38 Mean :18.02
## 3rd Qu.: 8.000 3rd Qu.: 9.937 3rd Qu.:18.14 3rd Qu.:18.69
## Max. :41.000 Max. :15.879 Max. :20.97 Max. :20.97
##
## con.dens_curt duration season
## Min. :0.000 Min. : 9.00 Length:2083
## 1st Qu.:0.000 1st Qu.:24.00 Class :character
## Median :1.260 Median :32.00 Mode :character
## Mean :1.114 Mean :27.43
## 3rd Qu.:2.000 3rd Qu.:35.00
## Max. :3.448 Max. :35.00
## ## tree
mod.snow.tree.grow1 <- lmer(grow.snd ~ site + treatment + A.con_curt + A.het_curt +
(1|quad.unique) + (1|census) + (1|sp.),
data=subset(snowdat.g_adult, growth.form == 'tree'))
mod.snow.tree.grow2 <- lmer(grow.snd ~ site + treatment*duration+
(1|quad.unique) + (1|sp.),
data=subset(snowdat.g_adult, growth.form == 'tree'))
mod.snow.tree.grow3 <- lmer(grow.snd ~ site + treatment*season+
(1|quad.unique) + (1|sp.) + (1|census),
data=subset(snowdat.g_adult, growth.form == 'tree'))
anova(mod.snow.tree.grow1, mod.snow.tree.grow2, mod.snow.tree.grow3)## refitting model(s) with ML (instead of REML)## Data: subset(snowdat.g_adult, growth.form == "tree")
## Models:
## mod.snow.tree.grow2: grow.snd ~ site + treatment * duration + (1 | quad.unique) +
## mod.snow.tree.grow2: (1 | sp.)
## mod.snow.tree.grow1: grow.snd ~ site + treatment + A.con_curt + A.het_curt + (1 |
## mod.snow.tree.grow1: quad.unique) + (1 | census) + (1 | sp.)
## mod.snow.tree.grow3: grow.snd ~ site + treatment * season + (1 | quad.unique) + (1 |
## mod.snow.tree.grow3: sp.) + (1 | census)
## Df AIC BIC logLik deviance Chisq Chi Df
## mod.snow.tree.grow2 9 3624.0 3671.7 -1803.0 3606.0
## mod.snow.tree.grow1 10 3433.3 3486.3 -1706.7 3413.3 192.7130 1
## mod.snow.tree.grow3 10 3432.7 3485.6 -1706.3 3412.7 0.6714 0
## Pr(>Chisq)
## mod.snow.tree.grow2
## mod.snow.tree.grow1 < 2.2e-16 ***
## mod.snow.tree.grow3 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod.snow.tree.grow.diags <- augment(mod.snow.tree.grow1)
qplot(data=mod.snow.tree.grow.diags, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
geom_hline(yintercept=0, linetype='dotted')## `geom_smooth()` using method = 'gam'
qplot(data = mod.snow.tree.grow.diags, x=treatment, y=.wtres, geom="boxplot") 
qqPlot(resid(mod.snow.tree.grow1))
summary(mod.snow.tree.grow1, correlation=FALSE)## Linear mixed model fit by REML ['lmerMod']
## Formula: grow.snd ~ site + treatment + A.con_curt + A.het_curt + (1 |
## quad.unique) + (1 | census) + (1 | sp.)
## Data: subset(snowdat.g_adult, growth.form == "tree")
##
## REML criterion at convergence: 3444.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8592 -0.5555 -0.0079 0.5554 3.3538
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.000000 0.0000
## sp. (Intercept) 0.004665 0.0683
## census (Intercept) 0.142553 0.3776
## Residual 0.588307 0.7670
## Number of obs: 1471, groups: quad.unique, 59; sp., 13; census, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -0.6298832 0.4244555 -1.484
## site1 0.1347980 0.0340059 3.964
## site2 -0.0168361 0.0313303 -0.537
## treatmentS 0.0436222 0.0417431 1.045
## A.con_curt 0.0005533 0.0072676 0.076
## A.het_curt 0.0333985 0.0211320 1.580# although treating site as a fixed effect make the model performing better, results do not change
# duration has no effect on seedling growth, there is no need to consider it in the following analyses.
## shrub
mod.snow.shrub.grow1 <- lmer(grow.snd ~ site + treatment + A.sum_curt+
(1|quad.unique) + (1|census) + (1|sp.),
data=subset(snowdat.g_adult, growth.form == 'shrub'))
mod.snow.shrub.grow2 <- lmer(grow.snd ~ site + treatment +
(1|quad.unique) + (1|census) + (1|sp.),
data=subset(snowdat.g_adult, growth.form == 'shrub'))
mod.snow.shrub.grow3 <- lmer(grow.snd ~ treatment*season +
(1|quad.unique) + (1|census) + (1|sp.),
data=subset(snowdat.g_adult, growth.form == 'shrub'))
anova(mod.snow.shrub.grow1, mod.snow.shrub.grow2, mod.snow.shrub.grow3)## refitting model(s) with ML (instead of REML)## Data: subset(snowdat.g_adult, growth.form == "shrub")
## Models:
## mod.snow.shrub.grow2: grow.snd ~ site + treatment + (1 | quad.unique) + (1 | census) +
## mod.snow.shrub.grow2: (1 | sp.)
## mod.snow.shrub.grow3: grow.snd ~ treatment * season + (1 | quad.unique) + (1 | census) +
## mod.snow.shrub.grow3: (1 | sp.)
## mod.snow.shrub.grow1: grow.snd ~ site + treatment + A.sum_curt + (1 | quad.unique) +
## mod.snow.shrub.grow1: (1 | census) + (1 | sp.)
## Df AIC BIC logLik deviance Chisq Chi Df
## mod.snow.shrub.grow2 8 1240.3 1275.6 -612.14 1224.3
## mod.snow.shrub.grow3 8 1239.5 1274.8 -611.76 1223.5 0.7582 0
## mod.snow.shrub.grow1 9 1237.5 1277.2 -609.74 1219.5 4.0374 1
## Pr(>Chisq)
## mod.snow.shrub.grow2
## mod.snow.shrub.grow3 <2e-16 ***
## mod.snow.shrub.grow1 0.0445 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod.snow.shrub.grow.diags <- augment(mod.snow.shrub.grow1)
qplot(data=mod.snow.shrub.grow.diags, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
geom_hline(yintercept=0, linetype='dotted')## `geom_smooth()` using method = 'loess'
qplot(data = mod.snow.shrub.grow.diags, x=treatment, y=.wtres, geom="boxplot") 
qqPlot(resid(mod.snow.shrub.grow1))
summary(mod.snow.shrub.grow1, correlation=FALSE)## Linear mixed model fit by REML ['lmerMod']
## Formula: grow.snd ~ site + treatment + A.sum_curt + (1 | quad.unique) +
## (1 | census) + (1 | sp.)
## Data: subset(snowdat.g_adult, growth.form == "shrub")
##
## REML criterion at convergence: 1240.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.7576 -0.4549 0.0047 0.4997 3.1070
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.00000 0.0000
## sp. (Intercept) 0.00000 0.0000
## census (Intercept) 0.05456 0.2336
## Residual 0.42203 0.6496
## Number of obs: 612, groups: quad.unique, 51; sp., 13; census, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.36225 0.60169 2.264
## site1 -0.01962 0.04158 -0.472
## site2 0.02836 0.03670 0.773
## treatmentS -0.02627 0.05697 -0.461
## A.sum_curt -0.07138 0.03263 -2.188##------------------------- individual species -------------------------------##
## And now pulling the species apart.
## First need to filter down to species that can be analysed - need enough
## reps (lets say 50) and variation in abundance (at least >10
sp.snowgrow.counts <- summarise(group_by(snowdat.g_adult, sp.),
n.indiv=length(grow.snd), n.quadrat=length(unique(grow.snd)))
sp.snowgrow.sel <- sp.snowgrow.counts$sp.[sp.snowgrow.counts$n.indiv > 50 & sp.snowgrow.counts$n.quadrat > 10]
sp.snowgrow.counts[sp.snowgrow.counts$sp. %in% sp.snowgrow.sel,]## # A tibble: 8 × 3
## sp. n.indiv n.quadrat
## <fctr> <int> <int>
## 1 ACPS 75 14
## 2 DEGL 68 24
## 3 EUAL 65 19
## 4 EUMA 96 26
## 5 EUVE 69 14
## 6 FRMA 793 24
## 7 PHSC 179 43
## 8 TIAM 379 13sp.snowgrow.sel2 <- c('ACPS','DEGL','EUAL','EUMA', 'EUVE','FRMA','PHSC','TIAM')
snowgrowmods <- sapply(sp.snowgrow.sel2, function(sp){
print(sp)
spdat <- filter(snowdat.g_adult, sp. == sp)
spdat <- droplevels(spdat)
mod <- lmer(grow.snd ~ site+ treatment + A.sum_curt + (1|quad.unique) + (1|census),
data=spdat)
}, simplify=FALSE)## [1] "ACPS"
## [1] "DEGL"
## [1] "EUAL"
## [1] "EUMA"
## [1] "EUVE"
## [1] "FRMA"
## [1] "PHSC"
## [1] "TIAM"## diagnostic plots
indmods.snowgrow.resids <- lapply(snowgrowmods, augment)
indmods.snowgrow.resids <- do.call('rbind', mapply(function(spnm, dat) {
dat$species <- spnm
return(dat)
}, spnm=as.list(names(indmods.snowgrow.resids)), dat=indmods.snowgrow.resids,
SIMPLIFY=FALSE))
ggplot(data=indmods.snowgrow.resids, aes(x=.fitted, y=.wtres, colour=treatment)) +
geom_jitter() + geom_smooth() + facet_wrap(~species, scales='free') +
geom_hline(yintercept=0, linetype='dashed')## `geom_smooth()` using method = 'loess'
library(lattice)
lapply(snowgrowmods, function(mod){
dotplot(ranef(mod, condVar=TRUE))})## $ACPS
## $ACPS$quad.unique
##
## $ACPS$census
##
##
## $DEGL
## $DEGL$quad.unique
##
## $DEGL$census
##
##
## $EUAL
## $EUAL$quad.unique
##
## $EUAL$census
##
##
## $EUMA
## $EUMA$quad.unique
##
## $EUMA$census
##
##
## $EUVE
## $EUVE$quad.unique
##
## $EUVE$census
##
##
## $FRMA
## $FRMA$quad.unique
##
## $FRMA$census
##
##
## $PHSC
## $PHSC$quad.unique
##
## $PHSC$census
##
##
## $TIAM
## $TIAM$quad.unique
##
## $TIAM$census
lapply(snowgrowmods, function(x) summary(x, correlation=FALSE))## $ACPS
## Linear mixed model fit by REML ['lmerMod']
## Formula: grow.snd ~ site + treatment + A.sum_curt + (1 | quad.unique) +
## (1 | census)
## Data: spdat
##
## REML criterion at convergence: 157
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3769 -0.5956 -0.1252 0.7792 2.0944
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 1.518e-16 1.232e-08
## census (Intercept) 2.306e-01 4.802e-01
## Residual 4.086e-01 6.392e-01
## Number of obs: 75, groups: quad.unique, 23; census, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.92470 2.19067 0.422
## site2 -0.26127 0.20462 -1.277
## site3 -0.46614 0.20549 -2.268
## treatmentS 0.10123 0.23310 0.434
## A.sum_curt -0.02566 0.12401 -0.207
##
## $DEGL
## Linear mixed model fit by REML ['lmerMod']
## Formula: grow.snd ~ site + treatment + A.sum_curt + (1 | quad.unique) +
## (1 | census)
## Data: spdat
##
## REML criterion at convergence: 141.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8994 -0.3190 -0.0513 0.4379 2.2436
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.0000 0.0000
## census (Intercept) 0.1721 0.4149
## Residual 0.3839 0.6196
## Number of obs: 68, groups: quad.unique, 10; census, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.86927 2.92275 1.324
## site2 0.32105 0.28425 1.129
## site3 -0.13271 0.39301 -0.338
## treatmentS 0.02295 0.17033 0.135
## A.sum_curt -0.21884 0.16525 -1.324
##
## $EUAL
## Linear mixed model fit by REML ['lmerMod']
## Formula: grow.snd ~ site + treatment + A.sum_curt + (1 | quad.unique) +
## (1 | census)
## Data: spdat
##
## REML criterion at convergence: 153.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.63460 -0.38641 -0.04121 0.64050 2.06026
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 1.717e-15 4.144e-08
## census (Intercept) 1.073e-01 3.276e-01
## Residual 5.709e-01 7.556e-01
## Number of obs: 65, groups: quad.unique, 9; census, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.01245 2.67990 1.124
## site2 -0.01935 0.27799 -0.070
## site3 -0.37808 0.38383 -0.985
## treatmentS -0.03388 0.79567 -0.043
## A.sum_curt -0.15730 0.15082 -1.043
##
## $EUMA
## Linear mixed model fit by REML ['lmerMod']
## Formula: grow.snd ~ site + treatment + A.sum_curt + (1 | quad.unique) +
## (1 | census)
## Data: spdat
##
## REML criterion at convergence: 190.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.22987 -0.45917 -0.03787 0.54778 2.72495
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 6.175e-18 2.485e-09
## census (Intercept) 7.848e-02 2.801e-01
## Residual 3.712e-01 6.093e-01
## Number of obs: 96, groups: quad.unique, 17; census, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 4.20258 1.67419 2.510
## site2 0.24947 0.17292 1.443
## site3 0.25304 0.28763 0.880
## treatmentS -0.11429 0.18321 -0.624
## A.sum_curt -0.23330 0.09523 -2.450
##
## $EUVE
## Linear mixed model fit by REML ['lmerMod']
## Formula: grow.snd ~ site + treatment + A.sum_curt + (1 | quad.unique) +
## (1 | census)
## Data: spdat
##
## REML criterion at convergence: 64.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.52899 -0.72941 0.00044 0.55533 2.95673
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.0000 0.0000
## census (Intercept) 0.3640 0.6033
## Residual 0.1102 0.3320
## Number of obs: 69, groups: quad.unique, 25; census, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -1.86794 1.52112 -1.228
## site2 -0.25097 0.13524 -1.856
## site3 -0.39067 0.13243 -2.950
## treatmentS 0.11792 0.13009 0.906
## A.sum_curt 0.09580 0.08444 1.135
##
## $FRMA
## Linear mixed model fit by REML ['lmerMod']
## Formula: grow.snd ~ site + treatment + A.sum_curt + (1 | quad.unique) +
## (1 | census)
## Data: spdat
##
## REML criterion at convergence: 1743
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9285 -0.6169 0.0257 0.5867 2.9782
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.0000 0.0000
## census (Intercept) 0.2182 0.4671
## Residual 0.5048 0.7105
## Number of obs: 793, groups: quad.unique, 52; census, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.20195 0.51213 0.394
## site2 -0.15881 0.06735 -2.358
## site3 -0.14253 0.07743 -1.841
## treatmentS 0.06777 0.05315 1.275
## A.sum_curt -0.01362 0.02699 -0.505
##
## $PHSC
## Linear mixed model fit by REML ['lmerMod']
## Formula: grow.snd ~ site + treatment + A.sum_curt + (1 | quad.unique) +
## (1 | census)
## Data: spdat
##
## REML criterion at convergence: 351.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.6747 -0.3637 -0.0267 0.4419 3.1665
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.0009391 0.03064
## census (Intercept) 0.0295244 0.17183
## Residual 0.3786474 0.61534
## Number of obs: 179, groups: quad.unique, 20; census, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.73684 0.90776 0.812
## site2 0.08695 0.10975 0.792
## site3 0.10830 0.14705 0.736
## treatmentS 0.04211 0.10331 0.408
## A.sum_curt -0.04240 0.05140 -0.825
##
## $TIAM
## Linear mixed model fit by REML ['lmerMod']
## Formula: grow.snd ~ site + treatment + A.sum_curt + (1 | quad.unique) +
## (1 | census)
## Data: spdat
##
## REML criterion at convergence: 1000.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9642 -0.6932 0.1655 0.5997 2.6766
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 8.748e-16 2.958e-08
## census (Intercept) 1.436e-02 1.198e-01
## Residual 7.914e-01 8.896e-01
## Number of obs: 379, groups: quad.unique, 55; census, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -0.52974 0.98108 -0.540
## site2 -0.15713 0.11993 -1.310
## site3 -0.37056 0.11569 -3.203
## treatmentS -0.01124 0.09531 -0.118
## A.sum_curt 0.03878 0.05464 0.710lapply(snowgrowmods, function(x) anova(x))## $ACPS
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 3.6090 1.80449 4.4167
## treatment 1 0.0614 0.06141 0.1503
## A.sum_curt 1 0.0175 0.01750 0.0428
##
## $DEGL
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 1.22278 0.61139 1.5925
## treatment 1 0.17727 0.17727 0.4617
## A.sum_curt 1 0.67330 0.67330 1.7538
##
## $EUAL
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 0.56452 0.28226 0.4944
## treatment 1 0.03197 0.03197 0.0560
## A.sum_curt 1 0.62099 0.62099 1.0878
##
## $EUMA
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 0.06733 0.03367 0.0907
## treatment 1 0.60827 0.60827 1.6385
## A.sum_curt 1 2.22815 2.22815 6.0019
##
## $EUVE
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 0.86165 0.43083 3.9078
## treatment 1 0.04710 0.04710 0.4272
## A.sum_curt 1 0.14190 0.14190 1.2871
##
## $FRMA
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 3.9208 1.96039 3.8832
## treatment 1 0.7907 0.79071 1.5663
## A.sum_curt 1 0.1286 0.12858 0.2547
##
## $PHSC
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 0.17711 0.088552 0.2339
## treatment 1 0.12116 0.121164 0.3200
## A.sum_curt 1 0.25772 0.257717 0.6806
##
## $TIAM
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## site 2 7.7853 3.8927 4.9186
## treatment 1 0.0702 0.0702 0.0888
## A.sum_curt 1 0.3985 0.3985 0.5036##------------------------ Other species ------------------------##
snowdat.g.other_adult <- snowdat.g_adult[!(snowdat.g_adult$sp. %in% sp.snowgrow.sel2),]
mod.snow.other.grow <- lmer(grow.snd ~ site + treatment+ A.sum_curt +
(1|quad.unique) + (1|census)+ (1|sp.), data=snowdat.g.other_adult)
mod.snow.other.grow.diags <- augment(mod.snow.other.grow)
qplot(data=mod.snow.other.grow.diags, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
geom_hline(yintercept=0, linetype='dotted')## `geom_smooth()` using method = 'loess'
qplot(data = mod.snow.other.grow.diags, x=treatment, y=.wtres, geom="boxplot") 
qqPlot(resid(mod.snow.other.grow))
summary(mod.snow.other.grow, correlation=FALSE)## Linear mixed model fit by REML ['lmerMod']
## Formula: grow.snd ~ site + treatment + A.sum_curt + (1 | quad.unique) +
## (1 | census) + (1 | sp.)
## Data: snowdat.g.other_adult
##
## REML criterion at convergence: 797.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.4727 -0.5532 0.0165 0.5332 3.0289
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.00232 0.04817
## sp. (Intercept) 0.00000 0.00000
## census (Intercept) 0.06168 0.24836
## Residual 0.50184 0.70840
## Number of obs: 359, groups: quad.unique, 45; sp., 18; census, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.40890 0.97006 0.422
## site1 0.01117 0.05755 0.194
## site2 -0.04173 0.05384 -0.775
## treatmentS -0.07238 0.08630 -0.839
## A.sum_curt -0.01542 0.05285 -0.292Snow removal did not affect the seedling growth.
Diversity analysis
library(vegan)## Loading required package: permute## This is vegan 2.4-4divers_woodsnow <- summarise(
group_by(snowdat.ag_adult, site, quadrat, quad.unique, treatment, census),
shannon=diversity(survs), simpson=diversity(survs, index='invsimpson'), species=specnumber(survs))
adult.ngbr.inf <- read.csv("data/adult.ngbr.inf.csv")
divers_woodsnow$A.sum <- adult.ngbr.inf$A.sum[match(divers_woodsnow$quad.unique, adult.ngbr.inf$quad.unique)]
summary(divers_woodsnow)## site quadrat quad.unique treatment census
## 1:111 Min. :102.0 1_1_102: 6 C:172 15sp:53
## 2:114 1st Qu.:203.0 1_1_104: 6 S:147 15fa:46
## 3: 94 Median :306.0 1_1_203: 6 16sp:43
## Mean :299.5 1_1_205: 6 16fa:59
## 3rd Qu.:404.0 1_1_206: 6 17sp:58
## Max. :508.0 1_1_210: 6 17fa:60
## (Other):283
## shannon simpson species A.sum
## Min. :0.0000 Min. :1.000 Min. :0.000 Min. :4487
## 1st Qu.:0.0000 1st Qu.:1.684 1st Qu.:1.000 1st Qu.:5321
## Median :0.8018 Median :2.286 Median :3.000 Median :5662
## Mean :0.8094 Mean : Inf Mean :3.028 Mean :5900
## 3rd Qu.:1.2730 3rd Qu.:3.341 3rd Qu.:4.000 3rd Qu.:6504
## Max. :2.1458 Max. : Inf Max. :9.000 Max. :9223
## divers_woodsnow$A.sum_curt <- divers_woodsnow$A.sum^(1/3)
# season
divers_woodsnow$season <- ifelse(divers_woodsnow$census=='15fa'|
divers_woodsnow$census=='16fa'|
divers_woodsnow$census=='17fa', 'fa', 'sp')
summary(divers_woodsnow)## site quadrat quad.unique treatment census
## 1:111 Min. :102.0 1_1_102: 6 C:172 15sp:53
## 2:114 1st Qu.:203.0 1_1_104: 6 S:147 15fa:46
## 3: 94 Median :306.0 1_1_203: 6 16sp:43
## Mean :299.5 1_1_205: 6 16fa:59
## 3rd Qu.:404.0 1_1_206: 6 17sp:58
## Max. :508.0 1_1_210: 6 17fa:60
## (Other):283
## shannon simpson species A.sum
## Min. :0.0000 Min. :1.000 Min. :0.000 Min. :4487
## 1st Qu.:0.0000 1st Qu.:1.684 1st Qu.:1.000 1st Qu.:5321
## Median :0.8018 Median :2.286 Median :3.000 Median :5662
## Mean :0.8094 Mean : Inf Mean :3.028 Mean :5900
## 3rd Qu.:1.2730 3rd Qu.:3.341 3rd Qu.:4.000 3rd Qu.:6504
## Max. :2.1458 Max. : Inf Max. :9.000 Max. :9223
##
## A.sum_curt season
## Min. :16.49 Length:319
## 1st Qu.:17.46 Class :character
## Median :17.82 Mode :character
## Mean :18.02
## 3rd Qu.:18.67
## Max. :20.97
## ggplot(divers_woodsnow, aes(x=treatment, y=shannon)) + geom_boxplot()
summarise(group_by(divers_woodsnow, treatment, census), meanH=mean(shannon),
seH=sd(shannon)/sqrt(length(shannon))) %>%
ggplot(aes(x=treatment, y=meanH, ymin=meanH-seH, ymax=meanH+seH, colour=census)) +
geom_point(position=position_dodge(width = .5))+
geom_errorbar(position=position_dodge(width = .5), width = 0.2) +
coord_trans(y="exp") +
labs(x="Snow removal", y="Effective number of species") 
###-------------------- Overall ---------------------------###
# shannon diversity
mod.wsnow.shan1 <- lmer(shannon ~ site + treatment + (1|quad.unique) + (1|census),
data=divers_woodsnow)
mod.wsnow.shan2 <- lmer(shannon ~ site +treatment +A.sum_curt +(1|quad.unique)+ (1|census),
data=divers_woodsnow)
mod.wsnow.shan3 <- lmer(shannon ~ treatment*season + (1|quad.unique) + (1|census),
data=divers_woodsnow)
anova(mod.wsnow.shan1, mod.wsnow.shan2, mod.wsnow.shan3)## refitting model(s) with ML (instead of REML)## Data: divers_woodsnow
## Models:
## mod.wsnow.shan1: shannon ~ site + treatment + (1 | quad.unique) + (1 | census)
## mod.wsnow.shan3: shannon ~ treatment * season + (1 | quad.unique) + (1 | census)
## mod.wsnow.shan2: shannon ~ site + treatment + A.sum_curt + (1 | quad.unique) +
## mod.wsnow.shan2: (1 | census)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mod.wsnow.shan1 7 290.37 316.72 -138.18 276.37
## mod.wsnow.shan3 7 287.63 313.99 -136.82 273.63 2.734 0 <2e-16
## mod.wsnow.shan2 8 291.63 321.76 -137.82 275.63 0.000 1 1
##
## mod.wsnow.shan1
## mod.wsnow.shan3 ***
## mod.wsnow.shan2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(mod.wsnow.shan3, type=c('p', 'smooth'))
qqPlot(resid(mod.wsnow.shan3))
sjp.lmer(mod.wsnow.shan3, type='fe')## Computing p-values via Kenward-Roger approximation. Use `p.kr = FALSE` if computation takes too long.
sjp.lmer(mod.wsnow.shan3, type='re', sort.est="sort.all" )## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.## Plotting random effects...## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.## Plotting random effects...
summary(mod.wsnow.shan3, correlation=FALSE)## Linear mixed model fit by REML ['lmerMod']
## Formula: shannon ~ treatment * season + (1 | quad.unique) + (1 | census)
## Data: divers_woodsnow
##
## REML criterion at convergence: 282.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.42405 -0.60104 -0.00628 0.59849 2.36740
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.18839 0.4340
## census (Intercept) 0.11358 0.3370
## Residual 0.08021 0.2832
## Number of obs: 319, groups: quad.unique, 60; census, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.97980 0.21230 4.615
## treatmentS -0.17670 0.12073 -1.464
## seasonsp -0.22362 0.27856 -0.803
## treatmentS:seasonsp -0.25068 0.06445 -3.889##----------------------- More diversity indices -------------------------##
### Invsimpson
range(divers_woodsnow$simpson)## [1] 1 Infdivers_woodsnow1 <- subset(divers_woodsnow, simpson != Inf)
mod.wsnow.simp <- lmer(simpson ~ site + treatment*season + (1|quad.unique)+(1|census),
data=divers_woodsnow1)
plot(mod.wsnow.simp, type=c('p', 'smooth'))
sjp.lmer(mod.wsnow.simp, type='fe')## Computing p-values via Kenward-Roger approximation. Use `p.kr = FALSE` if computation takes too long.
sjp.lmer(mod.wsnow.simp, type='re', sort.est="sort.all" )## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.## Plotting random effects...
qqPlot(resid(mod.wsnow.simp))
summary(mod.wsnow.simp, correlation=FALSE)## Linear mixed model fit by REML ['lmerMod']
## Formula: simpson ~ site + treatment * season + (1 | quad.unique) + (1 |
## census)
## Data: divers_woodsnow1
##
## REML criterion at convergence: 766.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.50704 -0.54750 -0.00479 0.49057 2.62465
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.6946 0.8334
## census (Intercept) 0.3738 0.6114
## Residual 0.4640 0.6812
## Number of obs: 300, groups: quad.unique, 59; census, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.7789 0.3914 7.099
## site1 0.4896 0.1633 2.998
## site2 0.1122 0.1633 0.687
## treatmentS -0.3765 0.2436 -1.545
## seasonsp -0.4360 0.5103 -0.854
## treatmentS:seasonsp -0.4229 0.1613 -2.621### Pielou's evenness (J)
table(divers_woodsnow$species)##
## 0 1 2 3 4 5 6 7 8 9
## 19 63 65 62 33 39 17 12 7 2divers_woodsnow2 <- subset(divers_woodsnow, species > 1)
# 20% of the data don't meet the criteria
divers_woodsnow2$evenness <- divers_woodsnow2$shannon/log(divers_woodsnow2$species)
mod.wsnow.even <- lmer(evenness ~ site + treatment*season + (1|quad.unique) + (1|census),
data=divers_woodsnow2)
plot(mod.wsnow.even, type=c('p', 'smooth'))
sjp.lmer(mod.wsnow.even, type='fe')## Computing p-values via Kenward-Roger approximation. Use `p.kr = FALSE` if computation takes too long.
sjp.lmer(mod.wsnow.even, type='re', sort.est="sort.all" )## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.## Plotting random effects...
qqPlot(resid(mod.wsnow.even))
summary(mod.wsnow.even, correlation=FALSE)## Linear mixed model fit by REML ['lmerMod']
## Formula: evenness ~ site + treatment * season + (1 | quad.unique) + (1 |
## census)
## Data: divers_woodsnow2
##
## REML criterion at convergence: -459.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.6649 -0.6106 0.1042 0.6724 2.4046
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.004242 0.06513
## census (Intercept) 0.003701 0.06083
## Residual 0.004832 0.06951
## Number of obs: 237, groups: quad.unique, 56; census, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.885054 0.038181 23.180
## site1 0.032898 0.013948 2.359
## site2 -0.051215 0.013906 -3.683
## treatmentS 0.005687 0.021367 0.266
## seasonsp 0.029418 0.051147 0.575
## treatmentS:seasonsp -0.002994 0.019567 -0.153###-------------------- Trees, shrubs -------------------------###
divers_snow_gf <- summarise(
group_by(snowdat.ag_adult, site, quadrat, quad.unique, growth.form, treatment, census),
shannon=diversity(survs), simpson=diversity(survs, index='invsimpson'), species=specnumber(survs))
divers_snow_gf$A.sum <- adult.ngbr.inf$A.sum[match(divers_snow_gf$quad.unique, adult.ngbr.inf$quad.unique)]
divers_snow_gf$A.sum_curt <- divers_snow_gf$A.sum^(1/3)
summary(divers_snow_gf)## site quadrat quad.unique growth.form treatment census
## 1:191 Min. :102.0 1_1_305: 12 shrub:247 C:279 15sp: 86
## 2:188 1st Qu.:203.0 1_1_306: 12 tree :266 S:234 15fa: 67
## 3:134 Median :306.0 1_1_308: 12 16sp: 56
## Mean :299.9 1_1_401: 12 16fa: 97
## 3rd Qu.:404.0 1_1_501: 12 17sp: 96
## Max. :508.0 2_1_102: 12 17fa:111
## (Other):441
## shannon simpson species A.sum
## Min. :0.0000 Min. :1.000 Min. :0.000 Min. :4487
## 1st Qu.:0.0000 1st Qu.:1.000 1st Qu.:1.000 1st Qu.:5335
## Median :0.5004 Median :1.800 Median :2.000 Median :5662
## Mean :0.4698 Mean : Inf Mean :1.883 Mean :5924
## 3rd Qu.:0.7595 3rd Qu.:2.667 3rd Qu.:3.000 3rd Qu.:6533
## Max. :1.7329 Max. : Inf Max. :6.000 Max. :9223
##
## A.sum_curt
## Min. :16.49
## 1st Qu.:17.47
## Median :17.82
## Mean :18.05
## 3rd Qu.:18.69
## Max. :20.97
## divers_snow_gf$season <- ifelse(divers_snow_gf$census=='15fa'|
divers_snow_gf$census=='16fa'|
divers_snow_gf$census=='17fa', 'fa', 'sp')
summary(divers_snow_gf)## site quadrat quad.unique growth.form treatment census
## 1:191 Min. :102.0 1_1_305: 12 shrub:247 C:279 15sp: 86
## 2:188 1st Qu.:203.0 1_1_306: 12 tree :266 S:234 15fa: 67
## 3:134 Median :306.0 1_1_308: 12 16sp: 56
## Mean :299.9 1_1_401: 12 16fa: 97
## 3rd Qu.:404.0 1_1_501: 12 17sp: 96
## Max. :508.0 2_1_102: 12 17fa:111
## (Other):441
## shannon simpson species A.sum
## Min. :0.0000 Min. :1.000 Min. :0.000 Min. :4487
## 1st Qu.:0.0000 1st Qu.:1.000 1st Qu.:1.000 1st Qu.:5335
## Median :0.5004 Median :1.800 Median :2.000 Median :5662
## Mean :0.4698 Mean : Inf Mean :1.883 Mean :5924
## 3rd Qu.:0.7595 3rd Qu.:2.667 3rd Qu.:3.000 3rd Qu.:6533
## Max. :1.7329 Max. : Inf Max. :6.000 Max. :9223
##
## A.sum_curt season
## Min. :16.49 Length:513
## 1st Qu.:17.47 Class :character
## Median :17.82 Mode :character
## Mean :18.05
## 3rd Qu.:18.69
## Max. :20.97
## ### Trees
# shannon diversity
mod.treesnow.shan <- lmer(shannon ~ site + treatment*season + (1|quad.unique) + (1|census),
data=subset(divers_snow_gf, growth.form=='tree'))
plot(mod.treesnow.shan, type=c('p', 'smooth'))
qqPlot(resid(mod.treesnow.shan))
sjp.lmer(mod.treesnow.shan, type='fe')## Computing p-values via Kenward-Roger approximation. Use `p.kr = FALSE` if computation takes too long.
sjp.lmer(mod.treesnow.shan, type='re', sort.est="sort.all" )## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.## Plotting random effects...
summary(mod.treesnow.shan, correlation=FALSE)## Linear mixed model fit by REML ['lmerMod']
## Formula: shannon ~ site + treatment * season + (1 | quad.unique) + (1 |
## census)
## Data: subset(divers_snow_gf, growth.form == "tree")
##
## REML criterion at convergence: 206.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.5386 -0.5750 0.0214 0.5469 2.6580
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.07709 0.2777
## census (Intercept) 0.10313 0.3211
## Residual 0.07707 0.2776
## Number of obs: 266, groups: quad.unique, 60; census, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.57533 0.19489 2.952
## site1 0.12769 0.05645 2.262
## site2 0.02275 0.05633 0.404
## treatmentS -0.01821 0.08564 -0.213
## seasonsp -0.29592 0.26661 -1.110
## treatmentS:seasonsp -0.14373 0.07048 -2.039## Invsimpson
divers_snow_gf1 <- subset(divers_snow_gf, simpson != Inf)
mod.treesnow.simp <- lmer(simpson ~ site + treatment*season + (1|quad.unique) + (1|census),
data=subset(divers_snow_gf1, growth.form=='tree'))
plot(mod.treesnow.simp, type=c('p', 'smooth'))
qqPlot(resid(mod.treesnow.simp))
sjp.lmer(mod.treesnow.simp, type='fe')## Computing p-values via Kenward-Roger approximation. Use `p.kr = FALSE` if computation takes too long.
sjp.lmer(mod.treesnow.simp, type='re', sort.est="sort.all" )## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.## Plotting random effects...
summary(mod.treesnow.simp, correlation=FALSE)## Linear mixed model fit by REML ['lmerMod']
## Formula: simpson ~ site + treatment * season + (1 | quad.unique) + (1 |
## census)
## Data: subset(divers_snow_gf1, growth.form == "tree")
##
## REML criterion at convergence: 462.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.92527 -0.55317 -0.08986 0.48331 3.11755
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.3505 0.5920
## census (Intercept) 0.1581 0.3976
## Residual 0.2660 0.5158
## Number of obs: 224, groups: quad.unique, 59; census, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.90319 0.26199 7.264
## site1 0.38336 0.12017 3.190
## site2 -0.11325 0.11998 -0.944
## treatmentS 0.01019 0.17933 0.057
## seasonsp -0.42786 0.33846 -1.264
## treatmentS:seasonsp -0.16242 0.15467 -1.050### Shrubs
# shannon diversity
mod.shrubsnow.shan <- lmer(shannon ~ site + treatment*season + (1|quad.unique) + (1|census),
data=subset(divers_snow_gf, growth.form=='shrub'))
plot(mod.shrubsnow.shan, type=c('p', 'smooth'))
qqPlot(resid(mod.shrubsnow.shan))
sjp.lmer(mod.shrubsnow.shan, type='fe')## Computing p-values via Kenward-Roger approximation. Use `p.kr = FALSE` if computation takes too long.
sjp.lmer(mod.shrubsnow.shan, type='re', sort.est="sort.all" )## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.## Plotting random effects...
summary(mod.shrubsnow.shan, correlation=FALSE)## Linear mixed model fit by REML ['lmerMod']
## Formula: shannon ~ site + treatment * season + (1 | quad.unique) + (1 |
## census)
## Data: subset(divers_snow_gf, growth.form == "shrub")
##
## REML criterion at convergence: 32.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.11544 -0.50524 0.01993 0.30473 3.13103
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.166842 0.40846
## census (Intercept) 0.006976 0.08352
## Residual 0.030990 0.17604
## Number of obs: 247, groups: quad.unique, 52; census, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.54115 0.09828 5.506
## site1 -0.02548 0.08027 -0.317
## site2 0.06135 0.08153 0.752
## treatmentS -0.36828 0.11839 -3.111
## seasonsp -0.03984 0.07472 -0.533
## treatmentS:seasonsp -0.07002 0.04584 -1.527## Invsimpson
mod.shrubsnow.simp <- lmer(simpson ~ site + treatment*season + (1|quad.unique) + (1|census),
data=subset(divers_snow_gf1, growth.form=='shrub'))
plot(mod.shrubsnow.simp, type=c('p', 'smooth'))
qqPlot(resid(mod.shrubsnow.simp))
sjp.lmer(mod.shrubsnow.simp, type='fe')## Computing p-values via Kenward-Roger approximation. Use `p.kr = FALSE` if computation takes too long.
sjp.lmer(mod.shrubsnow.simp, type='re', sort.est="sort.all" )## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.## Plotting random effects...
summary(mod.shrubsnow.simp, correlation=FALSE)## Linear mixed model fit by REML ['lmerMod']
## Formula: simpson ~ site + treatment * season + (1 | quad.unique) + (1 |
## census)
## Data: subset(divers_snow_gf1, growth.form == "shrub")
##
## REML criterion at convergence: 279.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9466 -0.3996 -0.0146 0.2641 3.8626
##
## Random effects:
## Groups Name Variance Std.Dev.
## quad.unique (Intercept) 0.55828 0.7472
## census (Intercept) 0.01539 0.1240
## Residual 0.08765 0.2961
## Number of obs: 237, groups: quad.unique, 51; census, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.88602 0.17157 10.993
## site1 -0.08559 0.14817 -0.578
## site2 0.12985 0.15021 0.864
## treatmentS -0.64812 0.21748 -2.980
## seasonsp -0.06676 0.11374 -0.587
## treatmentS:seasonsp -0.10162 0.07972 -1.275